setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM")
dir1="./bayes_pvalue_beta0/"
dir2="./bayes_BF/"

filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210316LIBD.eQTL处理/53K.add.GWAS.eQTL.DEG.motif.for.analysis.txt",header=T,sep="\t")
#file$id=paste(file$Chr,file$Start,sep=":")
file=file[!duplicated(file$id),]	#对TF去重
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")
i=1
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=file1[file1$unitID %in% as.character(filea$unitID),]

file1$normal_var_freq=as.numeric(gsub("%","",file1$normal_var_freq))/100
file1$tumor_var_freq=as.numeric(gsub("%","",file1$tumor_var_freq))/100

file1=data.frame(unitID=file1$unitID,normal_var_freq=file1$normal_var_freq,tumor_var_freq=file1$tumor_var_freq)

str=unlist(strsplit(group1[i],"_"))
names(file1)=c("unitID",paste(str[1],".VAF",sep = ""),paste(str[2],".VAF",sep = ""))
result=file1

for(i in 2:14){
fn1=paste0(dir1,group1[i],".bayes_p.txt")
file1=read.table(fn1,head=T,sep = "\t")
file1$unitID=paste(file1$chrom,file1$position,file1$ref,file1$var,sep=":")
file1=file1[file1$unitID %in% as.character(filea$unitID),]

file1$normal_var_freq=as.numeric(gsub("%","",file1$normal_var_freq))/100
file1$tumor_var_freq=as.numeric(gsub("%","",file1$tumor_var_freq))/100

file1=data.frame(unitID=file1$unitID,normal_var_freq=file1$normal_var_freq,tumor_var_freq=file1$tumor_var_freq)
str=unlist(strsplit(group1[i],"_"))
names(file1)=c("unitID",paste(str[1],".VAF",sep = ""),paste(str[2],".VAF",sep = ""))
result=merge(result,file1,by="unitID",all=T)
}

write.table(result,"20210321.H3K分析/117012.VAF.txt",quote=F,row.names = T,sep="\t")

affects=c("X1T","M7","M5","M1","M47","M49","M28","M27","M30","M29","M26","M25","M35","M36")
unaffects=c("X2B","M8","M6","M2","M48","M50","M18","M17","M20","M19","M22","M21","M40","M39")

sel1=paste0(affects,".VAF")
sel2=paste0(unaffects,".VAF")

result$affect.VAF=rowMeans(result[,sel1],na.rm = T)
result$unaffect.VAF=rowMeans(result[,sel2],na.rm = T)
result$chazhi=result$affect.VAF-result$unaffect.VAF
result$alt.group=ifelse(result$chazhi>0,"up","down")
rt=data.frame(unitID=result$unitID,result[,30:33])
write.table(result,"20210321.H3K分析/117012.vaf.up.down.20210321.txt",quote=F,row.names = T,sep="\t")
